#This file changes the order of the polynomial for the main results of Dahlgaard (2016)

load("agg_day_year_all.rdata")
load("n_day_year_all.rdata")

## commented out for replication files 

#  install.packages("xtable")
#  install.packages("ggplot2")
#  install.packages("rmeta")
#  install.packages("dplyr")
#  install.packages("rdrobust")

library(xtable)
library(ggplot2)
library(dplyr)
library(rmeta)
library(rdrobust)

data <- left_join(data_day_year, data_n_year,
                  by = c("days", "year"))

data_close <- 
  data %>%
  mutate(voted     = round(par_vote*n),
         abstained = round((1-par_vote)*n),
         treated   = days > 0) %>% 
  select(year, days, treated, voted, abstained)

data_voted <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["voted"]])) %>%
  mutate(voted = 1)

data_abstained <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["abstained"]])) %>%
  mutate(voted = 0)

data <- 
  rbind(data_voted, data_abstained) 


# array with coefficients for years and models 

coef_array <- 
  array(NA, dim = c(2, 6, 5)) 

for (i in 1:4){
  year_vec <- c(2009, 2013, 2014, 2015)[i]
  data_year <- 
    data %>%
    filter(year == year_vec ) 
  rd1  <- with(data_year, rdrobust(y = voted, x = days, p = 1))
  rd2  <- with(data_year, rdrobust(y = voted, x = days, p = 2))
  rd3  <- with(data_year, rdrobust(y = voted, x = days, p = 3))
  rd4  <- with(data_year, rdrobust(y = voted, x = days, p = 4))
  rd5  <- with(data_year, rdrobust(y = voted, x = days, p = 5))
  rd6  <- with(data_year, rdrobust(y = voted, x = days, p = 6))
  coef_array[ , 1, i] <- c(rd1$coef[1],  rd1$se[1])
  coef_array[ , 2, i] <- c(rd2$coef[1],  rd2$se[1])
  coef_array[ , 3, i] <- c(rd3$coef[1],  rd3$se[1])
  coef_array[ , 4, i] <- c(rd4$coef[1],  rd4$se[1])
  coef_array[ , 5, i] <- c(rd5$coef[1],  rd5$se[1])
  coef_array[ , 6, i] <- c(rd6$coef[1],  rd6$se[1])
  print(i)
}

# Pool coefficients using fixed effects meta analysis
for (i in 1:6){
  coef_array[, i , 5] <- c(meta.summaries(coef_array[ 1, i, 1:4],
                                          coef_array[ 2, i, 1:4])$summary,
                           meta.summaries(coef_array[ 1, i, 1:4],
                                          coef_array[ 2, i, 1:4])$se.summary)
}

# create plot with effect estimates in different years

coef_data <- data_frame(est = 100*rev(coef_array[1,,5]),
                        se  = 100*rev(coef_array[2,,5]),
                        place = 1:6)

lab_vals <- rev(c("1st",
                 "2nd",
                 "3rd",
                 "4th",
                 "5th",
                 "6th"))

coef_plots_lc <- 
  ggplot(data = coef_data, aes(x = est,
                        y = place, 
                        xmin = est + qnorm(0.025) * se,
                        xmax = est + qnorm(0.975) * se) ) +
  geom_errorbarh(height = 0 ) + 
  geom_point(size = 3) + 
  theme_classic() + 
  xlab("Effect size in percentage points") +
  scale_y_continuous(name = "Order of polynomial",
                     breaks = 1:6,
                     labels = lab_vals) +
  geom_vline(xintercept = 0, linetype = 3, alpha = 0.5) +
  theme(axis.text = element_text(size = 16),
        axis.title = element_text(size = 20)) 
last_plot()
